function [Y X] = DataGeneration(model,n,m,psi,parameter)
% Generate correlated fixed effects with mean zero and covar matrix omega
rho =-.25; mu = zeros(n,2); omega = [1 rho;rho, 1]; log_fixed_effects = mvnrnd(mu,omega); 
alpha = exp(log_fixed_effects(:,1)); 
gamma = exp(log_fixed_effects(:,2));
% generate covariates
% first (binary) covariate:
threshold = -norminv(sqrt(1/2))*(1-rho);
temp = log(alpha)-log(gamma)>sqrt(2)*threshold; x2 = (temp*ones(1,m)).*(ones(n,1)*temp'); 
% second (continuous) covariate
mu0 = 0; beta0 = 0; rho0 = beta0/sqrt(1+beta0^2); sigma0 = 1; mean0 = mu0+sigma0*rho0*sqrt(2/pi); mu0 = -sigma0*rho0*sqrt(2/pi)+1;
mu1 = 0; beta1 = 3; rho1 = beta1/sqrt(1+beta1^2); sigma1 = 1; mean1 = mu1+sigma1*rho1*sqrt(2/pi); mu1 = -sigma1*rho1*sqrt(2/pi)-1;
u0 = randn(n,m); v0 = randn(n,m); w0 = rho0*u0+sqrt(1-rho0^2)*v0; z0 = w0.*(u0>=0)-w0.*(u0<0); xx0 = mu0 + sigma0*z0;
u1 = randn(n,m); v1 = randn(n,m); w1 = rho1*u1+sqrt(1-rho1^2)*v1; z1 = w1.*(u1>=0)-w1.*(u1<0); xx1 = mu1 + sigma1*z1;
x1 = xx0.*(x2==0)+xx1.*(x2==1);
X{1} = x1; X{2} = x2;
% generate condition-mean function
varphi = exp(x1*psi(1)+x2*psi(2)).*(alpha*ones(1,m)).*(ones(n,1)*gamma');
% generate outcome variable
if strcmp(model,'Poisson')==1,
   Y = random('Poisson',varphi); % poisson with mean = var = varphi
end
if strcmp(model,'Negbin2')==1,
    Y = random('nbin',parameter,parameter./(parameter+varphi)); % negbin2  with mean = varphi and var = mean(1+parameter*varphi)
end
if strcmp(model,'Lnormal')==1,
    if parameter==1, variance = 1             ; end % log-normal errors are standard normal
    if parameter==2, variance =   1./varphi   ; end % Poisson-type variance
    if parameter==3, variance = 1+1./varphi   ; end % Negbin2-type variance with theta = 1;
    if parameter==4, variance =   1./varphi.^2; end % outcomes are homoskedastic outcomes
    
    sigma2 = log(1+variance); mu = -sigma2/2; e = mu+sqrt(sigma2).*randn(n,m); epsilon = exp(e); 
    Y = varphi.*epsilon; % log-normal model  
    
    Y = round(Y); % introduce rounding error
    
end
if strcmp(model,'Mixture')==1,
    M = random('nbin',parameter,parameter./(parameter+varphi)); Y = zeros(n,m);
    for i = 1:n,
        for j = 1:m,
            Z = chi2rnd(1,M(i,j),1); 
            Y(i,j) = sum(Z);
        end
    end
    
    
end





    